Combined density-functional and dynamical cluster quantum Monte Carlo 
calculations for three-band Hubbard models for hole-doped cuprate superconductors 
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Using a combined local density functional theory (DFT-LDA) and quantum Monte Carlo (QMC) 
dynamic cluster approximation approach, the parameter dependence of the superconducting transi- 
tion temperature T c of several single-layer hole-doped cuprate superconductors with experimentally 
very different T cmzx is investigated. The parameters of two different three-band Hubbard models 
are obtained using the LDA and the downfolding TVth-order muffin-tin orbital technique with N = 
and 1 respectively. QMC calculations on 4-site clusters show that the d-wave transition tempera- 
ture T c depends sensitively on the parameters. While the A=1MT0 basis set which reproduces all 
three pda bands leads to a d-wave transition, the iV=0 set which merely reproduces the LDA Fermi 
surface and velocities does not. 



I. INTRODUCTION 

Despite intense experimental and theoretical efforts, 
an understanding of high-temperature superconductiv- 
ity (HTSC) in hole-doped cuprate materials is elusive. 
While the materials are increasingly well characterized 1 -, 
a firm theoretical linking of the superconducting transi- 
tion temperature T c to details of the underlying atom- 
istic and electronic structures remains a grand challenge 
in condensed-matter theory. 

Recent advances in quantum cluster theories have given 
insight into the two-dimensional (2D) one-band Hubbard 
model, the most commonly adopted many-electron model 
for these materials. At low temperatures and appro- 
priate hole concentrations, the model develops the req- 
uisite strong d x 2_ y 2 order and superconducting ground 
state 2 -^ due to magnetically driven pairing£££ If this 
model captures sufficient physics to describe real mate- 
rials, the magnetic interactions and resultant T c should 
be moderated by details of the actual materials and their 
electronic structures. 

In parallel with the Hubbard investigations, the elec- 
tronic structure of HTSCs have been extensively stud- 
ied using density functional theory^ (DFT). Although 
structural properties are well reproduced, conventional 
DFT-local density approximation (LDA) calculations fail 
to describe the undoped insulating ground stated ' 10 ' 11 ' 12 . 
Nevertheless, DFT calculations agree on a universal elec- 
tronic structure in these materials: the low-energy elec- 
tronic degrees of freedom are primarily the pda anti- 
bonding Op and Gud x 2_ y 2 orbitals in the C11O2 layer, 
and these bands have been parameterize d 13 ' 14 i 15 ' 16 . For 
all optimally and over-doped materials, the Fermi sur- 
faces (FSs) measured by angle-resolved photoemission 
(ARPES)i£ agree surprisingly well with detailed LDA 
predictions. Even the existence and distinct in-plane 
dispersion of the splitting between the two FS sheets 



in bi-layered cupratesi 3 - and the fc z -dispcrsion in body- 
centered tetragonal single-layered materials^, have re- 
cently been experimentally confirme d 17 i 19 . Finally, the 
LDA conduction-band parameter r ~ — t'/i, 13 which 
gives the material-dependence of the FS shape and has 
the same origin as t± (k||) has been found to correlate 
positively with T c 
established. 
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but causal links have not been 



In this paper we combine LDA-DFT calculations of the 
cuprates with quantum cluster calculations of the tran- 
sition temperatures T c of the 2D Hubbard model. Wc 
study the three- rather than the one-band model because 
the most localized Cu d x 2_ y 2-\ikc orbital describing the 
LDA conduction band is so extended^ that using merely 
the on-site Coulomb repulsion in a one-band model is not 
justified. However, including also the O x p x , and O y p y 
orbitals in the basis set, localizes the Cu d x 2_ y 2 orbital 
to the extent that the corresponding three-band Hubbard 
Hamiltonian appears to be a valid model. DFT calcula- 
tions are used to consistently obtain the parameters for 
the five single-layer materials HgBa2Cu04, T^P^CuOg, 
TlBaLaCu0 5 , La 2 Cu0 4 , and Ca 2 Cu0 2 Cl2, for which 
T cmax = 90 K, 85 K, 52 K, 40 K, and 26 K,respectively. 
To determine T c we use the dynamic cluster approxima- 
tion (DCA) with a finite-temperature quantum Monte 
Carlo (QMC) cluster solve r 21 i 22 i 23 and the previously 
calculated DFT parameters 2 ^. In principle, these cal- 
culations use no experimental input and are therefore a 
stringent test of both the density functional and quan- 
tum cluster methods, as well as the form of the under- 
lying model. As the first study of this type we aim to 
address the following questions: (1) What is the mag- 
nitude of T c variation in the Hubbard model following 
the LDA+DCA scheme, and is this variation realistic? 
(2) Are there parameters found by LDA-DFT beyond 
those typically considered in Hubbard-like schemes that 
are particularly important for determining T c in these 
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materials**. 



II. DENSITY FUNCTIONAL CALCULATIONS 

We approximate the LDA potential for the stoichiomet- 
ric (undoped) cuprates by a superposition of spherically- 
symmetric, overlapping potential wells and then con- 
struct the basis set of three orbitals per cell by downfold- 
ing within multiple-scattering theory at the LDA Fermi 
energy, €f- Such an orbital is constructed to have the 
following properties: (1) It solves Schrodinger's differ- 
ential equation at throughout the solid, i.e. in all 
partial-wave channels, except for a kink at the muffin-tin 
(MT) spheres in the Curf^-^, O x p x , and O y p y chan- 
nels. (2) It has no Cud x 2_ y 2, O x p x , or O y p y character 
inside any MT sphere other than the one in which the 
orbital is centered and has its own character. Hence, the 
orbital is chosen to vanish (with a kink) in the channels of 
the other orbitals, and this makes it maximally localized. 
Pictures are presented in Ref.— . This basis set of kinked 
partial waves^ provides Bloch solutions of Schrodinger's 
equation with errors proportional to e (k) - ej?, and en- 
ergy bands with errors proportional to [e (k) — €f] 2 due 
to the variational principle. The LDA FS and veloci- 
ties are thus given correctly, as can clearly be seen from 
Fig-Ill For the solutions of Schrodinger's equation at ep, 
the kinks cancel out. 

Instead of using kinked partial waves at as basis 
functions, we could have constructed iVth-order muffin- 
tin orbitals (TVMTOs)^ which for a mesh of energies, 
eo, ■-, ejVj yield wavefunctions with errors proportional to 
[e (k) — eo] .. [e (k) — ejy]- In that way, the three bands 
can be made to reproduce the LDA bands over a wider 
energy range, specifically the range set by Udd ~ 10 eV. In 
addition to the above-mentioned N=0 set, we shall also 
consider an A=l set with the second energy, e\, chosen 
near the bottom of the pda bonding band, i.e. 7 — 8 eV 
below eg = ep. As seen in Fig. QJ this basis set accounts 
for the LDA pda bonding, non-bonding, and antibonding 
bands, at the same time as it reproduces the LDA FS and 
velocities. However, its orbitals are slightly less localized. 

Symmetrical orthormalization of the three AMTOs fi- 
nally yields the orbital representation in which H is ex- 
pressed. The on-site elements of Hlda are the orbital 
energies, Cq DA , an< ^ the off-site elements are the integrals 
for hopping, t^ 1 , between orbital a on site ij and orbital 
(3 on site Im in the same CuC>2 layer. Inter-layer hopping 
is neglected in this study. The site indices for the d or- 
bital are integers, since it is on a cubic lattice, while those 
for the x (y) orbital are |0 (0^) plus integers. The or- 
bitals, centered on the Cu and O sites, are orthonormal 
and real by construction. On-site energies and hopping 
integrals are, therefore, real and symmetric. The Fourier 
components of iJ(k)LDA arc, for example iJ^k) — 
e ( i + 2t^QQ(cosfc I + cosfc a )+4t^Jcosfc K cosfc a + .... Table 
U gives the short notation used for the hopping integrals. 

The values of the hopping integrals for the N = ba- 



sis set are shown in Fig. [2] They have reasonably short 
range and are dominated by the usual t p d and t pp . The 
values ~ 0.9 eV of t p d are considerably smaller than the 
conventional value 1.5 e V 13 ' 14 i 15 ' 16 describing the width 
~ 4y/2t p d ~ 8.5 eV of the pda anti-bonding and bonding 
bands. Whereas t p d and most other hopping integrals 
are seen to be fairly independent of the material (only 
the one with apical-Cl is a bit smaller), t pp is not; it in- 
creases with the observed T c max . This is the conduction- 
band trend found previously^ and explained as p x to p y 
hopping via a high-energy (e s ) , Cu-centered, axial hy- 
brid consisting of Cu4s, Cu 3d 3z 2_ 1 , apical oxygen 2p z 
and axial cation orbitals, all stacked perpendicular to 
the layer. If the energy of this axial orbital is increased 

t 2 

(e.g. by moving apical oxygen closer to Cu), t pp ~ - ^ 
decreases^ 3 -. Within that axial model, t' pp = t pp , but 
Fig. [5] shows that this is not true: t' pp vanishes for the 
two high-T c cuprates, and for the three low-T c cuprates 
the sign of t' pp is opposite to that of t pp . The main reason 
is that hopping via the in-layer (and therefore material- 
independent) Cu 4p x orbital contributes to t' pp , but not 
to t pp , and opposes the hopping via the axial orbital^. 
Also material-independent hopping via Oj, p x orbitals in- 
fluences t' pp , and causes a sizeable t pp . So t pp and t' pp ex- 
hibit the material's trend. Finally, tdd proceeds mainly 
via the diffuse O x 3d 3x 2_ 1 orbital, lying 50 eV above ep, 
and t' pd proceeds mainly via polarization of the cation. 
In summary, (1) diffuse high-energy (e 7 ) orbitals make 
sizeable contributions cx i;* / (ep — e 7 ) to the LDA hop- 
ping integrals and (2) the energy e s of the axial orbital, 
here downfolded into the tails of the oxygen orbitals, is 
the essential material-dependent parameter. 

The Cu on-site Coulomb energy is Udd ~ 9.5 eV in all 
five cuprates as found by constrained LDA?- calculations 
with the LMTO-ASA method. The radius of the Cu 
sphere is adjusted to 1.32 A, such that the d-character 
in the upper half of the LDA conduction band, the d 
hole-count hd° A , is the same as that obtained from the 
three-band Hlda ■ This ensures that the Cud a .2_ y 2 par- 
tial wave truncated outside the atomic sphere is simi- 
lar to the Cnd x 2_ y 2 partner of the three-orbital AMTO 
set. For the four cuprates with apical oxygen, we find: 
h]f A = 0.46 and for the diagonal elements of H^da- 
e d° A (hd) ~ e p ~ 0.45 eV. Since correlation effects are al- 
ready taken into account at the mean-field level in the 
LDA, we must include a double-counting correction pro- 
portional to the deviation, /i^ DA — \ , of the d hole-count 
from that (^) of the d 9 — > d 1Q transition-state. The cor- 
rected orbital-energy difference is then: A = e p — = 
e-d + Udd ~e p = e^ A (h L d DA ) - e p + (h L d DA - \)U dd « 
0.0 eV, where e refers to the hole and e to the electron 
representation. The commonly assumed value is, how- 
ever, A ~ 3 eV i 14 ' 15 ! 16 Previous constrained LDA cal- 
culations for La2Cu0 4 14 ' 15 gave 3eV because they used 
the LMTO-ASA total d electron-count of 9.24 to deduce: 
^!f DA = 0.76, and then found the double-counting cor- 
rection to be ~2.5eV. However, integrating to the top of 
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TABLE I: Relationship between site index and short-form notation used for hopping integrals 



the conduction band, we find not 10, but 9.70 d electrons, 
which is consistent with /i^ DA = 0.46i2! 

Unfortunately, current many-body treatments fail to 
reproduce the insulating behavior at half-filling, unless 
A exceeds 3eV, and this is commonly felt to be unac- 
ceptable. We therefore empirically set A = 3.25 eV in 
Hlda, i.e. we increased e^ DA — e P by 3.25 eV, but kept our 
hopping integrals unchanged^ This does not completely 
ruin the agreement between the experimental and LDA 
FS shapes, but it weakens the trend: For Tl 2 Ba 2 CuC>6 
the effective —t'/t is reduced from 0.33^ to 0.22, with 
the experimental^ value being 0.28, while for La 2 CuC>4 
the reduction is merely from 0.17 to 0.16. In all five 
cases, the A-shift causes a 20% reduction of the effective 
conduction bandwidth 8t. 

With the N=l basis set, which describes the three 
LDA pda bands over the energy range Udd, rather than 
merely the antibonding band near the LDA Fermi level 
(see Fig.s[T]and|4]), the values of the hopping integrals for 
HgBa 2 Cu0 4 (90 K) and La 2 Cu0 4 (40 K) are as shown in 
FigJ21 Now t p d is increased to values much closer to the 
conventional one a 13 i 14 ' 15 i 16 and t pp is increased to 0.90 
eV. By having to span a wider energy-range, the A^=l 
orbitals are somewhat less localized, and consequently 
have somewhat longer-ranged hoppings, than the N=0 
orbitals. This also masks the material's trend in individ- 
ual hopping integrals, although it is of course present 
in the shape of the antibonding band near ef- Now 

e d DA ~ e P = °- 67 and °- 95 eV for HgBa 2 Cu0 4 (90 K) 
and La 2 Cu0 4 , respectively, but for the reason mentioned 
above, we shall set A to 3.25 eV. 



III. DYNAMIC CLUSTER APPROXIMATION 
CALCULATIONS 

To solve the 3-band Hubbard Hamiltonian H we use 
the DC A 22 i 23 (for a review, see Ref^). In this method 
we map the lattice model onto a periodic cluster of size 
L c x L c , embedded into a self-consistently determined 
mean-field background. Correlations up to a range ~L C 
are treated explicitly while longer-ranged correlations are 
treated at a mean field level. We solve the cluster prob- 
lem using QMC^i, which does not introduce further sig- 
nificant approximations. Calculations on large clusters at 
low temperatures become prohibitively (exponentially) 
expensive due to the QMC Fermion sign problem. T c is 
determined via the diverging d-wave pair-field suscepti- 
bility obtained over a series of calculations at progres- 
sively lower temperatures. We check for earlier diver- 
gences in other angular momentum channels. 

Due to the large computational cost of a parametric 



study using QMC, we have performed calculations on 4- 
site clusters (L c =2) at 15% hole doping, which is near 
optimal in real materials. The 4-site cluster is the small- 
est for which a d-wave order parameter is allowed topo- 
logically, and corresponds to a mean-field result 2 -. Single- 
band calculations on 4-site clusters have shown that the 
phase diagram of these clusters shows general agreement 
with HTSC— . Converged calculations on clusters of up 
to 26 sites - for a single set of parameters - find that T c of 
4-site clusters is over-estimated by a factor ^2^ i.e. the 
small clusters exhibit larger pairing correlations. Hence, 
the presence of d-wave order in 4-site clusters at low tem- 
perature does not confirm the existence of such order in 
larger clusters, while the absence of d-wave order strongly 
indicates an absence of this order in the thermodynamic 
limit. In all of our DC A calculations we consistently used 
the calculated hoppings t and not, for example, the orig- 
inal LDA dispersion. 

To establish the existence of a variation in T c in the 
three-band model we performed an initial parametric 
study using only the nearest-neighbor hopping integrals 
tpd and tpp. Leaving one of these fixed at the value calcu- 
lated with the 7V=0 basis set for HgBa 2 Cu0 4 , our highest 
T cmax material, we varied the other. As shown in Fig|5l 
a d-wave transition was obtained over the entire range 
of parameters studied. Increasing either parameter in- 
creased T c , and dT c /dt p d ~ 6.3 x dT c /dt pp . The increase 
of T c with increased t p d may be understood in terms of 
changes to the fundamental energy scale. Within the 
range of studied materials, however, the calculated varia- 
tion of tpd is rather insignificant due to very similar Cu-0 
bond lengths while t pp varies systematically, with larger 
values corresponding to materials with larger T cmax , a 
trend reproduced by our QMC calculations. 

In Fig JS] we show the calculated inverse d pair- field sus- 
ceptibility as a function of temperature for HgBa 2 Cu0 4 
with the N=0 basis set when all hoppings are included, 
as well as subsets of hoppings. Compared to the inverse 
susceptibility of the t p d-t pp only calculation (A), which 
yields a moderate T c ~ 17meV, when including all the 
hoppings (B) the inverse susceptibility is reduced at high 
temperatures, but reduces less quickly at lower temper- 
atures. Therefore, any transition for the true LDA N~Q 
hoppings must occur at much lower temperatures than 
for the tpd-tpp only case. The d hole occupancy at low 
temperature is hd ~0.76 compared to ~0.78 for the tpd- 
tpp case. For the complete set (B) there appears to be 
no d-wave transition at moderate temperatures. Due 
to the increasing computational cost for lower temper- 
atures we cannot completely exclude the possibility of 
a very-low-temperature transition, but it is certain that 
any T c is significantly reduced from the simpler tpd-t pp 
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only case. When calculations are performed for all five 
materials (not shown), we also find no apparent d-wave 
transitions. These results demonstrate that, surprisingly, 
T c of the three-band Hubbard model is a strong function 
of the hopping parameters beyond the nearest-neighbors. 

To investigate the cause of the T c reduction we system- 
atically surveyed the effect of varying each hopping pa- 
rameter to extract dT c / dt^™ . In (C) we see that adding 
tdd and t' dd to t p d and t pp is not what suppresses T c , but 
adding t' pd does, as seen in (D). Although all hopping 
parameters modify T c , variation of t' pd changes T c most 
dramatically. This hopping proceeds mainly via polariza- 
tion of the cation and is — O.lOeV for all five materials. 
Changing the sign of t' pd , but keeping all other hoppings 
realistic, even produces a significant enhancement of T c 
as shown in (E). This artificial sign change profoundly 
changes the U=0 conduction band: the effective t de- 
creases by a factor 2 from the full LDA value, and —t'/t 
decreases from 0.34 to 0, i.e. this change is opposite to 
the emperical trencU^. 

We now repeat the calculations using the N=l basis set 
which reproduces all three LDA pda bands over a range 
of 10 eV ~ Udd (see Fig. s [4] and [3]) . The inverse d pair- 
field susceptibilities for HgBa2Cu04 are shown in Fig. [7] 
When all the hoppings are included (B), the susceptibility 
is increased at all temperatures compared to the N~0 
case, and in contrast to the previous results, there does 
appear to be a d-wave transition at very low temperature. 
This result clearly demonstrates a strong sensitivity of 
the many-body results on details of the treatment of the 
LDA data and Hubbard Hamiltonian. 

To further investigate the differences between the N=Q 
and N=l parameter sets, we performed N=l calcula- 
tions including only t pp and t p d- From FigJ7]we see that 
neglecting all hoppings except these two between nearest 
neighbors has little effect on the susceptibility (A). This 
result is in marked contrast with the one for the 7V=0 set, 
where the removal of long-ranged terms, particularly t' d , 
significantly increased T c . This is presumably connected 
with the fact that t' pd changes sign when going from N=Q 
to N=l. We also show the effect of decreasing t pp with 
t p d = 1.28 eV (C,D): in contrast to results for t p d = 0.89 
eV (Fig[5]) we find T c to increase. This does not contra- 
dict the emperical trend that T c max increases with t'/t 
because this does not translate into an increase with t pp 
for the N—l set. Comparisons of the inverse suscepti- 
bility for calculations performed with the full parameters 
sets for HgBa2CuC>4 (B) and La2CuC>4 (E) indicates that 
La2Cu04 will have the higher transition temperature, the 
reverse order compared to experiment. 

As an additional independent test on the choice of 
downfolding method we repeated the QMC calculations 
using HgBa2CuC>4 hopping parameters determined by 
a Wannier-function projection method^ which, like the 
N=l set, reproduced all three pda bands. Nevertheless, 
in this case we found no d-wave transition in the compu- 
tationally accessible temperature range (a low tempera- 
ture transition cannot be ruled out). 



The above results clearly demonstrate that the phase 
diagram of the three-band Hubbard model is quite sen- 
sitive to the choice of hopping integrals, even around the 
commonly accepted energy range of t p d ~ leV. Although 
a d-wave transition is found for the iV=lMTO basis set 
which reproduces all three pda LDA bands as well as 
the LDA FS and velocities, no transition is found for 
slightly different choices of downfolding approach, e.g. 
for Wannier-function projection of the three bands or for 
the A=0MTO basis which only reproduces the LDA FS 
and velocities. We have only investigated a single point 
on the phase diagram due to the computational expense 
and numerical difficulty of the current QMC and DCA 
techniques. Within the three-band Hamiltonian, refine- 
ment of the ill-determined A = e p — e<j is clearly required, 
as well as investigation of the effect of different hoppings 
on the spectral properties. It is also highly desirable to 
investigate larger clusters as well as more complex Hamil- 
tonians: for example, in LDA-DFT the d 3z 2_ r 2 band lies 
close to the Fermi energy in some of the HTSC mate- 
rials, suggesting that additional Cu degrees of freedom 
may be required. Unfortunately these investigations are 
currently precluded due to the computational cost and 
worsening Fermion sign-problem; wc hope they will be 
examined in future. 



IV. CONCLUSIONS 

In summary, we have obtained the parameters of three- 
band Hubbard models for a series of single-layer cupratc 
superconductors with varying T c max from downfolding 
either to the LDA conduction band or to all three 
pda bands. The transition temperature calculated us- 
ing DCA-QMC on 4-site clusters and increasing the 
small LDA value of e p — Cd to 3.25 eV is a moderate to 
very strong function of the hopping parameters. Even 
small hopping integrals beyond the nearest-neighbours 
can have a marked effect on the transition temperatures. 
These parameters are sensitive to the choice of down- 
folding technique and effective degree of Wannier local- 
ization. The present calculations yields superconductiv- 
ity for the NMTO basis set which reproduces all three 
pda LDA bands, but not for the one which reproduces 
merely the LDA Fermi surface and velocities. We hope 
that our findings will motivate future investigations and 
methodological development of more robust approaches 
for constructing and/or for solving realistic models of the 
cuprate superconductors. 
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FIG. 1: Comparison of the downfolded N=0 MTO 3-band (bold) and complete LDA bandstructures (dashed) of four undoped 
single-layer HTSC materials. The zero of energy is the Fermi level. 
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FIG. 2: (Color online) Calculated in- layer hoppings for five single-layer cuprates for N=0. The materials are ordered by T CI 
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FIG. 3: (Color online) Comparison of calculated in-layer hoppings for two single-layer cuprates for N=0 and N=l (see text) 
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FIG. 4: Comparison of the downfolded N=l 3-band (bold) and complete LDA bandstructures (dashed) of the single-layer 
HTSC materials HgBa2Cu04 and La2Cu04. The zero of energy is the Fermi level. 
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FIG. 5: (Color online) Calculated transition temperatures T c for the 3-band t p d-t pp model (see text). The temperature variation 
is shown for variations in t pp with t p d, held fixed (blue crosses), and vise versa (red plusses). The guide lines are a linear fit. 
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FIG. 6: (Color online) Inverse d-wave pair-field susceptibility (P d )as a function of temperature for different hopping-parameter 
sets calculated from HgBa2Cu04 with N=Q (see text). The lines are power law fits to the lowest 5 points in each series. 



4.5 





+ 


A- 


HgBa2Cu04 t pd and t pp 


4 


X 


B - 


HgBa2Cu04 All hoppings 


3.5 


* 


C- 


HgBa2Cu04 t pd and t pp =0.7eV 




□ 


D - 


HgBa2Cu04 t pd and t pp =0.5eV ...qj'., 


3 


■ 


E - 


La2Cu04 All hoppings .••;>'*"*!#* 



■d2.5 h 
Q_ 

^ 2 ■ 



1.5 ■ 

1 - 

0.5 ■ 






HgBa 2 Cu0 4 and La 2 Cu0 4 
N=1 Parameters 



0.02 0.04 0.06 

Temperature (eV) 



0.08 



0.1 



FIG. 7: (Color online) Inverse d-wave pair-field susceptibility (P~ )as a function of temperature for different hopping-parameter 
sets calculated from HgBa2Cu04 and for La2Cu04 with iV=l (see text). The lines are power law fits to the lowest 5 points in 
each series. 



